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The main goal of electronic structure methods is to solve the Schrodinger equation 
for the electrons in a molecule or solid, to evaluate the resulting total energies, forces, 
response functions and other quantities of interest. In this paper we describe the 
basic ideas behind the main electronic structure methods such as the pseudopotential 
and the augmented wave methods and provide selected pointers to contributions that 
are relevant for a beginner. We give particular emphasis to the Projector Augmented 
Wave (PAW) method developed by one of us, an electronic structure method for ab- 
initio molecular dynamics with full wavefunctions. We feel that it allows best to 
show the common conceptional basis of the most widespread electronic structure 
methods in materials science. 

I. INTRODUCTION 

The methods described below require as input only the charge and mass of the nuclei, the 
number of electrons and an initial atomic geometry. They predict binding energies accurate 
within a few tenths of an electron volt and bond-lengths in the 1-2 percent range. Currently, 
systems with few hundred atoms per unit cell can be handled. The dynamics of atoms can 
be studied up to tens of pico-seconds. Quantities related to energetics, the atomic structure 
and to the ground-state electronic structure can be extracted. 

In order to lay a common ground and to define some of the symbols, let us briefly touch 
upon the density functional theory^^i^S. It maps a description for interacting electrons, a 
nearly intractable problem, onto one of non-interacting electrons in an effective potential. 
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Within density functional theory, the total energy is written as 
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Here, |\l/„) are one-particle electron states, /„ are the state occupations, n(r) = 
(r)\l/„(r) is the electron density and Z{r) = —J^R^R^i'^ ~ ^r) is the nuclear 
charge density density expressed in electron charges. is the atomic number of a nucleus 
at position R/j. It is implicitly assumed that the infinite self- interaction of the nuclei is 
removed. The exchange and correlation functional contains all the difficulties of the many- 
electron problem. The main conclusion of the density functional theory is that E^^ is a 
functional of the density. 

We use Dirac's bra and ket notation. A wavefunction \E'„, corresponds to a ket I^Pn), the 
complex conjugate wave function corresponds to a bra and a scalar product 

J (i^r\E'* (r)^m(r) is written as Vectors in the 3-d coordinate space are indicated 

by boldfaced symbols. Note that we use R as position vector and R as atom index. 
In current implementations, the exchange and correlation functional E;j.c[n{r)] has the form 



where F^c is a parameterized function of the density and its gradients. Such functionals are 
called gradient corrected. In local spin density functional theory, Fxc furthermore depends 
on the spin density and its derivatives. A review of the earlier developments has been given 
by Pari^. 

The electronic ground state is determined by minimizing the total energy functional -E'[\I'n] 
of Eq. ^ at a fixed ionic geometry. The one-particle wavefunctions have to be orthogonal. 
This constraint is implemented with the method of Lagrange multipliers. We obtain the 
ground state wavefunctions from the extremum condition for 



with respect to the wavefunctions and the Lagrange multipliers Am,n- The extremum con- 
dition for the wavefunctions has the form 
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where H = —-^-V"^ + Vcs{t) is the effective one-particle Hamilton operator. The effective 



2m, 

potential depends itself on the electron density via 

where fixd^) = '^^|n(r)'"^^ functional derivative of the exchange and correlation func- 

tional. 

After a unitary transformation that diagonalizes the matrix of Lagrange multipliers Am,n, 
we obtain the Kohn-Sham equations. 

iJl^O = l^n)en (4) 

The one-particle energies e„ are the eigenvalues of A„ ^ 

The remaining one-electron Schrodinger equations, namely the Kohn-Sham equations given 
above, still pose substantial numerical difficulties: (1) in the atomic region near the nucleus, 
the kinetic energy of the electrons is large, resulting in rapid oscillations of the wavefunction 
that require fine grids for an accurate numerical representation. On the other hand, the 
large kinetic energy makes the Schrodinger equation stiff, so that a change of the chemical 
environment has little effect on the shape of the wavefunction. Therefore, the wavefunction 
in the atomic region can be represented well already by a small basis set. (2) In the bonding 
region between the atoms the situation is opposite. The kinetic energy is small and the 
wavefunction is smooth. However, the wavefunction is flexible and responds strongly to the 
environment. This requires large and nearly complete basis sets. 

Combining these different requirements is non-trivial and various strategies have been de- 
veloped. 

• The atomic point of view has been most appealing to quantum chemists. Basis func- 
tions that resemble atomic orbitals are chosen. They exploit that the wavefunction in 
the atomic region can be described by a few basis functions, while the chemical bond 
is described by the overlapping tails of these atomic orbitals. Most techniques in this 
class are a compromise of, on the one hand, a well adapted basis set, where the basis 
functions are difficult to handle, and on the other hand numerically convenient basis 
functions such as Gaussians, where the inadequacies are compensated by larger basis 
sets. 
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• Pseudopotentials regard an atom as a perturbation of the free electron gas. The 
most natural basis functions are plane-waves. Plane waves basis sets are in principle 
complete and suitable for sufficiently smooth wavef unctions. The disadvantage of the 
comparably large basis sets required is offset by their extreme numerical simplicity. Fi- 
nite plane-wave expansions are, however, absolutely inadequate to describe the strong 
oscillations of the wavefunctions near the nucleus. In the pseudopotential approach 
the Pauli repulsion of the core electrons is therefore described by an effective potential 
that expels the valence electrons from the core region. The resulting wavefunctions 
are smooth and can be represented well by plane-waves. The price to pay is that all 
information on the charge density and wavefunctions near the nucleus is lost. 

• Augmented wave methods compose their basis functions from atom-like wavefunctions 
in the atomic regions and a set of functions, called envelope functions, appropriate 
for the bonding in between. Space is divided accordingly into atom-centered spheres, 
defining the atomic regions, and an interstitial region in between. The partial solutions 
of the different regions, are matched at the interface between atomic and interstitial 
regions. 

The projector augmented wave method is an extension of augmented wave methods and the 
pseudopotential approach, which combines their traditions into a unified electronic structure 
method. 

After describing the underlying ideas of the various approaches let us briefiy review the 
history of augmented wave methods and the pseudopotential approach. We do not discuss 
the atomic-orbital based methods, because our focus is the PAW method and its ancestors. 

II. AUGMENTED WAVE METHODS 

The augmented wave methods have been introduced in 1937 by Slatei»^ and were later 
modified by Korringa^, Kohn and Rostokker— . They approached the electronic structure 
as a scattered-electron problem. Consider an electron beam, represented by a plane-wave, 
traveling through a solid. It undergoes multiple scattering at the atoms. If for some energy, 
the outgoing scattered waves interfere destructively, a bound state has been determined. 
This approach can be translated into a basis set method with energy and potential dependent 
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basis functions. In order to make the scattered wave problem tractable, a model potential 
had to be chosen: The so-called muffin-tin potential approximates the true potential by a 
constant in the interstitial region and by a spherically symmetric potential in the atomic 
region. 

Augmented wave methods reached adulthood in the 1970s: O.K. Andersen^ showed that 
the energy dependent basis set of Slater's APW method can be mapped onto one with 
energy independent basis functions, by linearizing the partial waves for the atomic regions 
in energy. In the original APW approach, one had to determine the zeros of the determinant 
of an energy dependent matrix, a nearly intractable numerical problem for complex systems. 
With the new energy independent basis functions, however, the problem is reduced to the 
much simpler generalized eigenvalue problem, which can be solved using efficient numerical 
techniques. Furthermore, the introduction of well defined basis sets paved the way for full- 
potential calculation32^. In that case the muffin-tin approximation is used solely to define 
the basis set while the matrix elements {xi\H\xj) of the Hamiltonian are evaluated 
with the full potential. 

In the augmented wave methods one constructs the basis set for the atomic region by solving 
the radial Schrodinger equation for the spheridized effective potential 



e,r = 



_2me 

as function of energy. Note that a partial wave 0^,m(e, r) is an angular momentum eigenstate 
and can be expressed as a product of a radial function and a spherical harmonic. The energy 
dependent partial wave is expanded in a Taylor expansion about some reference energy e^^£ 

0^,m(e, r) = (j)u,e,m{r) + (e - e^^i)4>^^i^rn{r) + 0((e - e,,,^)^) 



where 0i//,m(r) = <Pi,mi^u,i,i^)- The energy derivative of the partial wave 0y(r) — ^^^^ 



de 



solves the equation 

-^2 



(r) = (f)u,i,m{r) 



2me 

Next one starts from a regular basis set, such as plane- waves, Gaussians or Hankel func- 
tions. These basis functions are called envelope functions \xi)- Within the atomic region 
they are replaced by the partial waves and their energy derivatives, such that the resulting 
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wavefunction is continuous and differentiable. 

_R R,e,rn 

Ofi{r) is a step function that is unity within the augmentation sphere centered at R/j and 
zero elsewhere. The augmentation sphere is atom-centered and has a radius about equal to 
the covalent radius. This radius is called the muffin-tin radius, if the spheres of neighboring 
atoms touch. These basis functions describe only the valence states; the core states are 
localized within the augmentation sphere and are obtained directly by radial integration of 
the Schrodinger equation within the augmentation sphere. 

The coefficients aR^e^rn,i and bR^£^m,i are obtained for each \xi) as follows: The envelope 
function is decomposed around each atomic site into spherical harmonics multiplied by 
radial functions. 

Xiir) = ^MR/,m,i(|r - R/j|)Y£,„(r - Rr) (6) 

Analytical expansions for plane-waves, Hankel functions or Gaussians exist. The radial parts 
of the partial waves (f)u,R,i,m and 4>u,R,e,m are matched with value and derivative to UR/,m,i(|r|), 
which yields the expansion coefficients aR^£^rn,i and bii/^rn,i- 

If the envelope functions are plane- waves, the resulting method is called the linear augmented 
plane-wave (LAPW) method. If the envelope functions are Hankel functions, the method is 
called linear muffin-tin orbital (LMTO) method . 

A good review of the LAPW method^ has been given by Singh4^. Let us now briefly men- 
tion the major developments of the LAPW method: SolerSfl introduced the idea of additive 
augmentation: While augmented plane-waves are discontinuous at the surface of the aug- 
mentation sphere if the expansion in spherical harmonics in Eq.Elis truncated, Soler replaced 
the second term in Eq. Elby an expansion of the plane- wave with the same angular momen- 
tum truncation as in the third term. This dramatically improved the convergence of the 
angular momentum expansion. Singb^ introduced so-called local orbitals, which are non- 
zero only within a muffin-tin sphere, where they are superpositions of (p and (p functions from 
different expansion energies. Local orbitals substantially increase the energy transferability. 
Sjostedl)^ relaxed the condition that the basis functions are differentiable at the sphere ra- 
dius. In addition she introduced local orbitals, which are confined inside the sphere, and 
that also have a kink at the sphere boundary. Due to the large energy-cost of kinks, they 
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will cancel, once the total energy is minimized. The increased variational degree of freedom 
in the basis leads to a dramatically improved plane-wave convergence^. 
The second variant of the linear methods is the LMTO method^. A good introduction into 
the LMTO method is the book by Skriver^^. The LMTO method uses Hankel functions 
as envelope functions. The atomic spheres approximation (ASA) provides a particularly 
simple and efficient approach to the electronic structure of very large systems. In the ASA 
the augmentation spheres are blown up so that their volume are equal to the total volume 
and the first two terms in Eq. Elare ignored. The main deficiency of the LMTO-ASA method 
is the limitation to structures that can be converted into a closed packed arrangement of 
atomic and empty spheres. Furthermore energy differences due to structural distortions are 
often qualitatively incorrect. Full potential versions of the LMTO method, that avoid these 
deficiencies of the ASA have been developed. The construction of tight binding orbitals as 
superposition of muffin-tin orbitals^ showed the underlying principles of the empirical tight- 
binding method and prepared the ground for electronic structure methods that scale linearly 
instead of with the third power of the number of atoms. The third generation LMTO^ allows 
to construct true minimal basis sets, which require only one orbital per electron-pair for 
insulators. In addition they can be made arbitrarily accurate in the valence band region, so 
that a matrix diagonalization becomes unnecessary. The first steps towards a full-potential 
implementation, that promises a good accuracy, while maintaining the simplicity of the of the 
LMTO-ASA method are currently under way. Through the minimal basis-set construction 
the LMTO method offers unrivaled tools for the analysis of the electronic structure and has 
been extensively used in hybrid methods combining density functional theory with model 
Hamiltonians for materials with strong electron correlationsti^ 



III. PSEUDOPOTENTIALS 

Pseudopotentials have been introduced to (1) avoid describing the core electrons explicitely 
and (2) to avoid the rapid oscillations of the wavefunction near the nucleus, which normally 
require either complicated or large basis sets. 

The pseudopotential approach traces back to 1940 when C. Herring invented the orthog- 
onalized plane- wave method^°. Later, Phillips-^ and Antoncik^ replaced the orthogonality 
condition by an effective potential, which mimics the Pauli-repulsion by the core electrons 
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and thus compensates the electrostatic attraction by the nucleus. In practice, the potential 
was modified, for example, by cutting off the singular potential of the nucleus at a certain 
value. This was done with a few parameters that have been adjusted to reproduce the 
measured electronic band structure of the corresponding solid. 

Hamann, Schliiter and Chiang^*^^ showed in 1979 how pseudopotentials can be constructed 
in such a way, that their scattering properties are identical to that of an atom to first 
order in energy. These first-principles pseudopotentials relieved the calculations from the 
restrictions of empirical parameters. Highly accurate calculations have become possible 
especially for semiconductors and simple metals. An alternative approach towards first- 
principles pseudopotentials^ preceeded the one mentioned above. 



A. The idea behind Pseudopotential construction 

In order to construct a first-principles pseudopotential, one starts out with an all-electron 
density-functional calculation for a spherical atom. Such calculations can be performed 
efficiently on radial grids. They yield the atomic potential and wavefunctions (f)e^rn{T^)- Due 
to the spherical symmetry, the radial parts of the wavefunctions for different magnetic 
quantum numbers m are identical. 

For the valence wavefunctions one constructs pseudo wavefunctions \4>e,m)'- There are nu- 
merous way3^'22»2^i^ to construct the pseudo wavefunctions. They must be identical to 
the true wave functions outside the augmentation region, which is called core-region in the 
context of the pseudopotential approach. Inside the augmentation region the pseudo wave- 
function should be node-less and have the same norm as the true wavefunctions, that is 
{4>e,m\4>e,m) = {(pe,m\(pe,m) (compare Figure C}. 

From the pseudo wavefunction, a potential ue{r) can be reconstructed by inverting the 
respective Schrodinger equation. 

This potential M^(r) (compare Figure QJ, which is also spherically symmetric, differs from 



V + ue{r) - ee.^ 



1 If ~ 
(r) = ^ u,{r) = e + ^— - ■ — -VVAn^(r) 



one main angular momentum i to the other. 

Next we define an effective pseudo Hamiltonian 

^2 'PS, ^ /■ ,n ,n(r') + Z(r') /r / m x 

= -2^.^ + 4^ y |r-r1 + "-""I--)!' ■■) 
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FIG. 1: Illustration of the pseudopotential concept at the example of the 3s wavefunction of Si. 
The solid line shows the radial part of the pseudo wavefunction The dashed line corresponds 

to the all-electron wavefunction <^i^ra which exhibits strong oscillations at small radii. The angular 
momentum dependent pseudopotential (dash-dotted line) deviates from the all-electron one Wg// 
(dotted line) inside the augmentation region. The data are generated by the fhi98PP code^. 

and determine the pseudopotentials such that the pseudo Hamiltonian produces the 
pseudo wavefunctions, that is 

<(r) = u,{r) - ^ / ^^'^^ + l^''^ - l^.m^)l r) (7) 
This process is called "unscreening" . 

Z(r) mimics the charge density of the nucleus and the core electrons. It is usually an atom- 
centered, spherical Gaussian that is normalized to the charge of nucleus and core of that 
atom. In the pseudopotential approach, Z^ij) it does not change with the potential. The 
pseudo density n(r) = ^^fn^^rSj)^ rSj) is constructed from the pseudo wavefunctions. 
In this way we obtain a different potential for each angular momentum channel. In order to 
apply these potentials to a given wavefunction, the wavefunction must first be decomposed 
into angular momenta. Then each component is applied to the pseudopotential for the 
corresponding angular momentum. 
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The pseudopotential defined in this way can be expressed in a semi-local form 



Si\v 










) 



The local potential v{r) only acts on those angular momentum components, not included in 
the expansion of the pseudopotential construction. Typically it is chosen to cancel the most 
expensive nonlocal terms, the one corresponding to the highest physically relevant angular 
momentum. 

The pseudopotential is non-local as it depends on two position arguments, r and r'. The 
expectation values are evaluated as a double integral 



The semi-local form of the pseudopotential given in Eq. |S1 is computationally expensive. 
Therefore, in practice one uses a separable form of the pseudopotential^i^SiS. 

-1 _ 



V 



Thus the projection onto spherical harmonics used in the semi-local form of Eq. IHlis replaced 
by a projection onto angular momentum dependent functions \v^''^(j)i). 

The indices i and j are composite indices containing the atomic-site index R, the angular 
momentum quantum numbers i,m and an additional index a. The index a distinguishes 
partial waves with otherwise identical indices R, i, m, as more than one partial wave per site 
and angular momentum is allowed. The partial waves may be constructed as eigenstates to 
the pseudopotential v^'^ for a set of energies. 

One can show that the identity of Eq. IHl holds by applying a wavefunction = J2i \4>i)ci 
to both sides. If the set of pseudo partial waves in Eq. IHlis complete, the identity is 
exact. The advantage of the separable form is that {(pv'^'^l is treated as one function, so that 
expectation values are reduced to combinations of simple scalar products 

B. The Pseudopotential total energy 

The total energy of the pseudopotential method can be written in the form 
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The constant Eseif is adjusted such that the total energy of the atom is the same for an 
all-electron calculation and the pseudopotential calculation. 

For the atom, from which it has been constructed, this construction guarantees that the 
pseudopotential method produces the correct one-particle energies for the valence states 
and that the wave functions have the desired shape. 

While pseudopotentials have proven to be accurate for a large variety of systems, there is 
no strict guarantee that they produce the same results as an all-electron calculation, if they 
are used in a molecule or solid. The error sources can be divided into two classes: 

• Energy transferability problems: Even for the potential of the reference atom, the 
scattering properties are accurate only in given energy window. 

• Charge transferability problems: In a molecule or crystal, the potential differs from 
that of the isolated atom. The pseudopotential, however, is strictly valid only for the 
isolated atom. 

The plane-wave basis set for the pseudo wavefunctions is defined by the shortest wave length 
A = 27r/|G| via the so-called plane wave cutoff Epw = ~2m^' often specified in 

Rydberg (1 Ry=i H^13.6 eV). The plane -wave cutoff is the highest kinetic energy of all 
basis functions. The basis-set convergence can systematically be controlled by increasing 
the plane-wave cutoff. 

The charge transferability is substantially improved by including a nonlinear core 
correction-^ into the exchange-correlation term of Eq. ^| Hamanni^ showed, how to con- 
struct pseudopotentials also from unbound wavefunctions. Vanderbillj^^^ generalized the 
pseudopotential method to non-normconserving pseudopotentials, so-called ultra-soft pseu- 
dopotentials, which dramatically improves the basis-set convergence. The formulation of 
ultra-soft pseudopotentials has already many similarities with the projector augmented wave 
method. Truncated separable pseudopotentials suffer sometimes from so-called ghost states. 
These are unphysical core-like states, which render the pseudopotential useless. These prob- 
lems have been discussed by GonzeiS,. Quantities such as hyperfine parameters that depend 
on the full wavefunctions near the nucleus, can be extracted approximately^. A good review 
about pseudopotential methodology has been written by Payne^ and Singb^. 
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In 1985 R. Car and M. Parrinello published the ab-initio molecular dynamics methodi*. 
Simulations of the atomic motion have become possible on the basis of state-of-the-art 
electronic structure methods. Besides making dynamical phenomena and finite tempera- 
ture effects accessible to electronic structure calculations, the ab-initio molecular dynamics 
method also introduced a radically new way of thinking into electronic structure methods. 
Diagonalization of a Hamilton matrix has been replaced by classical equations of motion for 
the wavefunction coefficients. If one applies friction, the system is quenched to the ground 
state. Without friction truly dynamical simulations of the atomic structure are performed. 
Using thermostatsi^'i^'^^i^ simulations at constant temperature can be performed. The Car- 
Parrinello method treats electronic wavefunctions and atomic positions on an equal footing. 

IV. PROJECTOR AUGMENTED WAVE METHOD 

The Car-Parrinello method had been implemented first for the pseudopotential approach. 
There seemed to be unsurmountable barriers against combining the new technique with 
augmented wave methods. The main problem was related to the potential dependent basis 
set used in augmented wave methods: the Car-Parrinello method requires a well defined and 
unique total energy functional of atomic positions and basis set coefficients. Furthermore 
the analytic evaluation of the first partial derivatives of the total energy with respect to wave 
functions, and atomic position, the forces, must be possible. Therefore, it was one of 

the main goals of the PAW method to introduce energy and potential independent basis sets 
that were as accurate as the previously used augmented basis sets. Other requirements have 
been: (1) The method should at least match the efficiency of the pseudopotential approach 
for Car-Parrinello simulations. (2) It should become an exact theory when converged and 
(3) its convergence should be easily controlled. We believe that these criteria have been met, 
which explains why the PAW method becomes increasingly wide spread today. 

A. Transformation theory 

At the root of the PAW method lies a transformation, that maps the true wavefunctions 
with their complete nodal structure onto auxiliary wavefunctions, that are numerically con- 
venient. We aim for smooth auxiliary wavefunctions, which have a rapidly convergent plane- 
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wave expansion. With such a transformation we can expand the auxihary wave functions 
into a convenient basis set such as plane-waves, and evaluate all physical properties after 
reconstructing the related physical (true) wavef unctions. 

Let us denote the physical one-particle wavefunctions as |\E'„) and the auxiliary wavefunctions 
as l^^n)- Note that the tilde refers to the representation of smooth auxiliary wavefunctions 
and n is the label for a one-particle state and contains a band index, a /c-point and a spin 
index. The transformation from the auxiliary to the physical wave functions is denoted by 

r. 

I^n)=r|^„) (11) 

Now we express the constrained density functional F of Eq. El in terms of our auxiliary 
wavefunctions 

F[r^„, A^,„] = E[T^^] - J2[{^n\r^r\^m) - 6n,m]An.,n (12) 

n,m 

The variational principle with respect to the auxiliary wavefunctions yields 

r^HT\^^) = rtr|^„)e„. (13) 

Again we obtain a Schrodinger-like equation (see derivation of Eq. 0]), but now the Hamilton 
operator has a different form, H = T^HT, an overlap operator O = T'^T occurs, and the 
resulting auxiliary wavefunctions are smooth. 

When we evaluate physical quantities we need to evaluate expectation values of an operator 
A, which can be expressed in terms of either the true or the auxiliary wavefunctions. 

(A) = J2 /n(^n|A|vl>„) = J2 fn{^n\T^AT\^^) (14) 
n n 

In the representation of auxiliary wavefunctions we need to use transformed operators A = 
T'^AT. As it is, this equation only holds for the valence electrons. The core electrons are 
treated differently as will be shown below. 

The transformation takes us conceptionally from the world of pseudopotentials to that of 
augmented wave methods, which deal with the full wavefunctions. We will see that our 
auxiliary wavefunctions, which are simply the plane-wave parts of the full wavefunctions, 
translate into the wavefunctions of the pseudopotential approach. In the PAW method the 
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auxiliary wavefunctions are used to construct the true wavefunctions and the total energy 
functional is evaluated from the latter. Thus it provides the missing link between augmented 
wave methods and the pseudopotential method, which can be derived as a well-defined 
approximation of the PAW method. 

In the original paper^, the auxiliary wavefunctions have been termed pseudo wavefunctions 
and the true wavefunctions have been termed all-electron wavefunctions, in order to make 
the connection more evident. We avoid this notation here, because it resulted in confusion 
in cases, where the correspondence is not clear-cut. 



Sofar, we have described how we can determine the auxiliary wave functions of the ground 
state and how to obtain physical information from them. What is missing, is a definition of 
the transformation operator T. 

The operator T has to modify the smooth auxiliary wave function in each atomic region, so 
that the resulting wavefunction has the correct nodal structure. Therefore, it makes sense 
to write the transformation as identity plus a sum of atomic contributions Sr 



R 

For every atom, Sr adds the difference between the true and the auxiliary wavefunction. 
The local terms Sr are defined in terms of solutions of the Schrodinger equation for 
the isolated atoms. This set of partial waves \(f)i) will serve as a basis set so that, near the 
nucleus, all relevant valence wavefunctions can be expressed as superposition of the partial 
waves with yet unknown coefficients. 



With i G -R we indicate those partial waves that belong to site R. 

Since the core wavefunctions do not spread out into the neighboring atoms, we will treat them 
differently. Currently we use the frozen-core approximation, which imports the density and 
the energy of the core electrons from the corresponding isolated atoms. The transformation 
T shall produce only wavefunctions orthogonal to the core electrons, while the core electrons 
are treated separately. Therefore, the set of atomic partial waves includes only valence 
states that are orthogonal to the core wavefunctions of the atom. 



B. Transformation operator 




(15) 




(16) 
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For each of the partial waves we choose an auxihary partial wave |</)j). The identity 

\cl>i) = (1 + 5^)10,) for leR 
Sr\~^.) = 10.) - \4>i) (17) 

defines the local contribution Sr to the transformation operator. Since 1 + Sr shall change 
the wavefunction only locally, we require that the partial waves and their auxiliary 
counter parts are pairwise identical beyond a certain radius Tc^r. 

(/)i(r) = 0i(r) for i e R and |r - Kr\ > r^^R (18) 

Note that the partial waves are not necessarily bound states and are therefore not normaliz- 
able, unless we truncate them beyond a certain radius r^R. The PAW method is formulated 
such that the final results do not depend on the location where the partial waves are trun- 
cated, as long as this is not done too close to the nucleus and identical for auxiliary and 
all-electron partial waves. 

In order to be able to apply the transformation operator to an arbitrary auxiliary wave- 
function, we need to be able to expand the auxiliary wavefunction locally into the auxiliary 
partial waves. 

^{r) = ^^i{r)ci = ^4)i{v){Pi\^!) for |r - R^| < r^.R (19) 

i&R ieR 

which defines the projector functions \pi). The projector functions probe the local character 
of the auxiliary wave function in the atomic region. Examples of projector functions are 
shown in Figure El From Eq. ^Jwe can derive Xligi? \'Pi){Pi\ = 1; which is valid within r^ij. 
It can be shown by insertion, that the identity Eq. holds for any auxiliary wavefunction 
|\E') that can be expanded locally into auxiliary partial waves if 

{Pi\(i)j)=kj for i.jeR (20) 

Note that neither the projector functions nor the partial waves need to be orthogonal among 
themselves. The projector functions are fully determined with the above conditions and a 
closure relation, which is related to the unscreening of the pseudopotentials (see Eq. 90 in^). 
By combining Eq. El and Eq. we can apply Sr to any auxiliary wavefunction. 

5^1^) = = - (21) 

i&R i&R 
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FIG. 2: Projector functions of the chlorine atom. Top: two s-type projector functions, middle: 
p-type, bottom: d-type. 

Hence the transformation operator is 

T=l + (22) 

i 

where the sum runs over all partial waves of all atoms. The true wave function can be 
expressed as 

1^) = 1^) + - = 1^) + - 1^^)) (23) 

i R 

with 

I^^r) = El'^^)(^^l^) (24) 

i€R 

\%) = J2\^^){pm (25) 

ieR 

In Fig. El the decomposition of Eq. |2S1 is shown for the example of the bonding p-cr state of 
the CI2 molecule. 

To understand the expression Eq. ESI for the true wave function, let us concentrate on 
different regions in space. (1) Far from the atoms, the partial waves are, according to Eq. IH^ 
pairwise identical so that the auxiliary wavefunction is identical to the true wavefunction, 
that is ^I/(r) = \E'(r). (2) Close to an atom R, however, the auxiliary wavefunction is, 
according to Eq. UHl identical to its one-center expansion, that is \^(r) = ^^)j(r). Hence 
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FIG. 3: Bonding p-a orbital of the CI2 molecule and its decomposition of the wavefunction into 
auxiliary wavefunction and the two one-center expansions. Top-left: True and auxiliary wave 
function; top-right: auxiliary wavefunction and its partial wave expansion; bottom-left: the two 
partial wave expansions; bottom-right: true wavefunction and its partial wave expansion. 

the true wavefunction ^^'(r) is identical to \I/^(r), which is built up from partial waves that 
contain the proper nodal structure. 

In practice, the partial wave expansions are truncated. Therefore, the identity of Eq.imidoes 
not hold strictly. As a result the plane-waves also contribute to the true wavefunction inside 
the atomic region. This has the advantage that the missing terms in a truncated partial wave 
expansion are partly accounted for by plane-waves, which explains the rapid convergence of 
the partial wave expansions. This idea is related to the additive augmentation of the LAPW 
method of Solei^. 

Frequently, the question comes up, whether the transformation Eq. |221of the auxiliary wave- 
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functions indeed provides the true wavefunction. The transformation should be considered 
merely as a change of representation analogous to a coordinate transform. If the total energy 
functional is transformed consistently, its minimum will yield auxiliary wavefunctions that 
produce the correct wave functions |\1'). 



C. Expectation values 

Expectation values can be obtained either from the reconstructed true wavefunctions or 
directly from the auxiliary wave functions 



n n=l 
= J2 fn{^n\r^Ar\^^) + (26) 



n=l 



where /„ are the occupations of the valence states and Nc is the number of core states. The 
first sum runs over the valence states, and second over the core states 
Now we can decompose the matrix element for a wavefunction \1/ into its individual contri- 
butions according to Eq. ESI 

R R! 

= {^\m+Y,{{'^R\A\^R)-{^R\A\^R) 
R 



part 1 

+ 5^((vl/^ - ^],\A\^ - ^],) + - ^'n\A\^'n ' ^r) 



part 2 

+ J2^^R-^r\^\^R'-^R') (27) 



R^R' 



part 3 

Only the first part of Eq. 123 is evaluated explicitly, while the second and third parts of 
Eq.EZIare neglected, because they vanish for sufficiently local operators as long as the partial 
wave expansion is converged: The function — vanishes per construction beyond its 
augmentation region, because the partial waves are pairwise identical beyond that region. 
The function \& — vanishes inside its augmentation region, if the partial wave expansion 
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is sufficiently converged. In no region of space both functions — and \^ — are 
simultaneously nonzero. Similarly the functions — from different sites are never non- 
zero in the same region in space. Hence, the second and third parts of Eq. |2ZI vanish for 
operators such as the kinetic energy ^^V^ and the real space projection operator |r)(r|, 
which produces the electron density. For truly nonlocal operators the parts 2 and 3 of Eq.lTTI 
would have to be considered explicitly. 

The expression, Eq. 1211 for the expectation value can therefore be written with the help of 
Eq.lZHas 



Nc 



n=l 



n 

Nc ^ 

= J2fn{^n\A\^n)+J2(^n\Am 

n n=l 

Nc,R 

R ijGR neR 

Nc,R 

R i,jeR neR 

where Di j is the one-center density matrix defined as 

D^,J = J2fn{^n\pJ){p^\^n) = J^iP^l^ n) fn{^ n\Pj) (29) 
n n 

The auxiliary core states, allow to incorporate the tails of the core wavefunction into 
the plane-wave part, and therefore assure, that the integrations of partial wave contributions 
cancel strictly beyond Vc- They are identical to the true core states in the tails, but are a 
smooth continuation inside the atomic sphere. It is not required that the auxiliary wave 
functions are normalized. 

Following this scheme, the electron density is given by 

n{r) = n{r) + J] (^(r) - n)j(r)) (30) 

R 

n(r) = E/n^:(r)^n(r)+ne 

n 

i,j&R 

: E A>Kr)0«(r)+^c,fl (31) 

i,j&R 



nifr 
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where is the core density of the corresponding atom and nc^R is the auxihary core density 
that is identical to outside the atomic region, but smooth inside. 

Before we continue, let us discuss a special point: The matrix element of a general operator 
with the auxiliary wavefunctions may be slowly converging with the plane-wave expansion, 
because the operator A may not be well behaved. An example for such an operator is the 
singular electrostatic potential of a nucleus. This problem can be alleviated by adding an 
"intelligent zero" : If an operator B is purely localized within an atomic region, we can use 
the identity between the auxiliary wavefunction and its own partial wave expansion 

= {^n\B\^n) - {K\B\^l) (32) 

Now we choose an operator B so that it cancels the problematic behavior of the operator 
A, but is localized in a single atomic region. By adding B to the plane-wave part and the 
matrix elements with its one-center expansions, the plane- wave convergence can be improved 
without affecting the converged result. A term of this type, namely v will be introduced in 
the next section to cancel the Coulomb singularity of the potential at the nucleus. 



D. Total Energy 

Like wavefunctions and expectation values also the total energy can be divided into three 
parts. 

E[^^,Rn] = E + J2{e],-E},) (33) 

R 

The plane-wave part E involves only smooth functions and is evahiated on equi-spaced grids 
in real and reciprocal space. This part is computationally most demanding, and is similar 
to the expressions in the pseudopotential approach. 



1 



^3^ f^s^, [n{r) + Z{r)][n{r') + Z{r')] 



2 Attcq J J |r - r' 

+ J dhv{r)fi{T) + E^e[n(r)] (34) 

Z{r) is an angular-momentum dependent core-like density that will be described in detail 
below. The remaining parts can be evaluated on radial grids in a spherical harmonics expan- 



21 



sion. The nodal structure of the wavefunctions can be properly described on a logarithmic 
radial grid that becomes very fine near nucleus, 

i,j€R ^ neR ^ 

+ 1.^ /^s, f,,^ M\r) + Z{r)][n\r') + Z{r')] 



2 47reo 

+ E,,[n\r)] (35) 

e'r = 



2me 



i,j€R 

1 

2 ' ivreo 



+ J d^rv{r)h\r) + E^^[n\r)] (36) 

The compensation charge density Z{r) = ^/?(r) is given as a sum of angular momentum 
dependent Gauss functions, which have an analytical plane-wave expansion. A similar term 
occurs also in the pseudopotential approach. In contrast to the norm-conserving pseudopo- 
tential approach, however, the compensation charge of an atom Z^ is non-spherical and 
constantly adapts to the instantaneous environment. It is constructed such that 



Zn{r) - n],ir) - Z^(r) (37) 



has vanishing electrostatic multi-pole moments for each atomic site. With this choice, the 
electrostatic potentials of the augmentation densities vanish outside their spheres. This is 
the reason that there is no electrostatic interaction of the one-center parts between different 
sites. 

The compensation charge density as given here is still localized within the atomic regions. A 
technique similar to an Ewald summation, however, allows to replace it by a very extended 
charge density. Thus we can achieve, that the plane-wave convergence of the total energy is 
not affected by the auxiliary density. 

The potential v = J2r ^R^ which occurs in Eqs. EH and EHl enters the total energy in the 
form of "intelligent zeros" described in Eq. IHS 

= J2fn [{^n\vR\^n) ' {^1,\vr\K)) = /n (^n | ^n) " ^^AkvR\4>j) (38) 

n n i,j&R 

The main reason for introducing this potential is to cancel the Coulomb singularity of 
the potential in the plane- wave part. The potential v allows to influence the plane- wave 
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convergence beneficially, without changing the converged result, v must be localized within 
the augmentation region, where Eq. holds. 

E. Approximations 

Once the total energy functional provided in the previous section has been defined, ev- 
erything else follows: Forces are partial derivatives with respect to atomic positions. The 
potential is the derivative of the non-kinetic energy contributions to the total energy with 
respect to the density, and the auxiliary Hamiltonian follows from derivatives with re- 

spect to auxiliary wave functions. The fictitious Lagrangian approach of Car and Parrinello^ 
does not allow any freedom in the way these derivatives are obtained. Anything else than 
analytic derivatives will violate energy conservation in a dynamical simulation. Since the 
expressions are straightforward, even though rather involved, we will not discuss them here. 
All approximations are incorporated already in the total energy functional of the PAW 
method. What are those approximations? 

• Firstly we use the frozen-core approximation. In principle this approximation can be 
overcome. 

• The plane-wave expansion for the auxiliary wavefunctions must be complete. The 
plane-wave expansion is controlled easily by increasing the plane-wave cutoff defined 
as Epw = I^^G^Q^. Typically we use a plane-wave cutoff of 30 Ry. 

• The partial wave expansions must be converged. Typically we use one or two partial 
waves per angular momentum {£, m) and site. It should be noted that the partial wave 
expansion is not variational, because it changes the total energy functional and not 
the basis set for the auxiliary wavefunctions. 

We do not discuss here numerical approximations such as the choice of the radial grid, since 
those are easily controlled. 

F. Relation to the Pseudopotentials 

We mentioned earlier that the pseudopotential approach can be derived as a well defined 
approximation from the PAW method: The augmentation part of the total energy AE = 



23 



E-^ — for one atom is a functional of the one-center density matrix Di j^ji defined in 
Eq. 1221 The pseudopotential approach can be recovered if we truncate a Taylor expansion 
of AE about the atomic density matrix after the linear term. The term linear to Dij is the 
energy related to the nonlocal pseudopotential. 

f)/\ TP 

AE{D,,) = AE(Dg) + $^(A,-A-5)^ + 0(A,-Al)' 

= E,,if + ^ f„,{^,,\vn^n) - I d^rv{v)n{v) + 0(A,, - A'P' (39) 

n ^ 

which can directly be compared to the total energy expression Eq. [TUlof the pseudopotential 
method. The local potential 'i;(r) of the pseudopotential approach is identical to the corre- 
sponding potential of the projector augmented wave method. The remaining contributions 
in the PAW total energy, namely E, differ from the corresponding terms in Eq. E3only in 
two features: our auxiliary density also contains an auxiliary core density, reflecting the non- 
linear core correction of the pseudopotential approach, and the compensation density Z{r) 
is non-spherical and depends on the wave function. Thus we can look at the PAW method 
also as a pseudopotential method with a pseudopotential that adapts to the instantaneous 
electronic environment. In the PAW method, the explicit nonlinear dependence of the total 
energy on the one-center density matrix is properly taken into account. 
What are the main advantages of the PAW method compared to the pseudopotential ap- 
proach? 

Firstly all errors can be systematically controlled so that there are no transferability errors. 
As shown by Watson^^ and Kresse^i, most pseudopotentials fail for high spin atoms such 
as Cr. While it is probably true that pseudopotentials can be constructed that cope even 
with this situation, a failure can not be known beforehand, so that some empiricism remains 
in practice: A pseudopotential constructed from an isolated atom is not guaranteed to be 
accurate for a molecule. In contrast, the converged results of the PAW method do not 
depend on a reference system such as an isolated atom, because PAW uses the full density 
and potential. 

Like other all-electron methods, the PAW method provides access to the full charge and spin 
density, which is relevant, for example, for hyperfine parameters. Hyperfine parameters are 
sensitive probes of the electron density near the nucleus. In many situations they are the 
only information available that allows to deduce atomic structure and chemical environment 
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of an atom from experiment. 

The plane-wave convergence is more rapid than in norm-conserving pseudopotentials and 
should in principle be equivalent to that of ultra-soft pseudopotentials^. Compared to the 
ultra-soft pseudopotentials, however, the PAW method has the advantage that the total 
energy expression is less complex and can therefore be expected to be more efficient. 
The construction of pseudopotentials requires to determine a number of parameters. As 
they influence the results, their choice is critical. Also the PAW methods provides some 
flexibility in the choice of auxiliary partial waves. However, this choice does not influence 
the converged results. 

G. Recent Developments 

Since the flrst implementation of the PAW method in the CP-PAW code, a number of 
groups have adopted the PAW method. The second implementation was done by the group 
of Holzwartb^^. The resulting PWPAW code is freely available^. This code is also used as 
a basis for the PAW implementation in the Ablnit project. An independent PAW code has 
been developed by Valiev and Weare^. Recently the PAW method has been implemented 
into the VASP code^^. The PAW method has also been implemented by W. Kromen into 
the ESTCoMPP code of Bliigel and Schroder. 

Another branch of methods uses the reconstruction of the PAW method, without taking 
into account the full wavefunctions in the energy minimization. Following chemist notation 
this approach could be termed "post-pseudopotential PAW" . This development began with 
the evaluation for hyperflne parameters from a pseudopotential calculation using the PAW 
reconstruction operator^"^ and is now used in the pseudopotential approach to calculate 
properties that require the correct wavefunctions such as hyperflne parameters. 
The implementation by Kresse and Jouberlj^ has been particularly useful as they had an 
implementation of PAW in the same code as the ultra-soft pseudopotentials, so that they 
could critically compare the two approaches with each other. Their conclusion is that both 
methods compare well in most cases, but they found that magnetic energies are seriously 
- by a factor two - in error in the pseudopotential approach, while the results of the PAW 
method were in line with other all-electron calculations using the linear augmented plane- 
wave method. As a short note, Kresse and Joubert incorrectly claim that their implemen- 
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tation is superior as it includes a term that is analogous to the non-linear core correction of 
pseudopotentials^^: this term however is already included in the original version in the form 
of the pseudized core density. 

Several extensions of the PAW have been done in the recent years: For applications in 
chemistry truly isolated systems are often of great interest. As any plane- wave based 
method introduces periodic images, the electrostatic interaction between these images can 
cause serious errors. The problem has been solved by mapping the charge density onto a 
point charge model, so that the electrostatic interaction could be subtracted out in a self- 
consistent mannev^. In order to include the influence of the environment, the latter was 
simulated by simpler force fields using the molecular-mechanics-quantum-mechanics (QM- 
MM) approach^. 

In order to overcome the limitations of the density functional theory several extensions have 
been performed. Bengone^ implemented the LDA+U approach into the CP-PAW code. 
Soon after this, Arnaud^ accomplished the implementation of the GW approximation into 
the CP-PAW code. The VASP-version of PAW^^ and the CP-PAW code have now been 
extended to include a non-coUinear description of the magnetic moments. In a non-collinear 
description the Schrodinger equation is replaced by the Pauli equation with two-component 
spinor wavefunctions 

The PAW method has proven useful to evaluate electric field gradient^^ and magnetic hy- 
perfine parameters with high accuracy^^. Invaluable will be the prediction of NMR chemical 
shifts using the GIPAW method of Pickard and Mauri^, which is based on their earlier 
worb^S. While the GIPAW is implemented in a post-pseudopotential manner, the extension 
to a self-consistent PAW calculation should be straightforward. An post-pseudopotential ap- 
proach has also been used to evaluate core level spectra^^ and momentum matrix elements^^. 
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